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Two specialized algorithms for the numerical integration of the equations of motion of a Brownian 
walker obeying detailed balance are introduced. The algorithms become symplectic in the appro- 
priate limits, and reproduce the equilibrium distributions to some higher order in the integration 
time step. Comparisons with other existing integration schemes are carried out both for static and 
dynamical quantities. 
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I. INTRODUCTION 

Stochastic processes are well known to be at the heart 
of many physical systems Q. Several approaches have 
hence been developed to understand the dynamics which 
is realized given specific models: among others, one of the 
most used ones is Monte-Carlo integration of the equa- 
tions of motion. 

The literature on the numerical integration of stochas- 
tic differential equations is huge: we will limit here to 
mention a couple of classical citations widely used in the 
physics community 0, 0, 0. Additional comments and 
references can be found in [f| . The numerical algorithms 
presented in these works are general and can be applied 
to basically any flow; however, they might not be the op- 
timal ones for cases when additional information about 
the details of the system under study are available. 

An important class for which dedicated algorithms can 
be derived is given by the equations of motion 



x 
v 



-jv + F(x) + f (t) 



(1) 



where is a random Gaussian noise, with zero average 
and standard deviation 

= WS(t- a). 

In the following, we will also use V[x), defined as F(x) = 
—V'(x). Note that although we are dealing here with 
only one Brownian walker, the algorithm we are going to 
show can be easily extended to the case when x and v are 
vectors, and Fi(x), the force acting on the i-th walker, is 
a function of all other walkers. 

The above equation is commonly found in the liquid 
state literature (for numerical schemes appropriate in the 
integration of the Brownian dynamics of a liquid, see 
among others [E ll> El EH ) an d some algorithms have 
been proposed, over the years, for its numerical integra- 
tion. 

To date, perhaps the most widely used algorithms for 
the integration of Eq. ^are the ones derived in @, where 
two algorithms have been proposed (see also references 
therein) : we will benchmark against one of them, and to 
this end, we will briefly review them below. Note that 



the system of Eq. ^ become s sy mplectic when 7^0 
and, until some recent works [TH [12. Il3| . this symplec- 
tic nature was not really exploited in deriving numerical 
schemes. 

The approach we will follow is to derive numerical al- 
gorithms having in mind two requirements: (i) the al- 
gorithm should become symplectic in the deterministic 
(T = 0) and frictionless (7 = 0) limit; and (ii) the nu- 
merical algorithm should reproduce as closely as possible 
the equilibrium distribution, when it is defined, of the 
system given by Eq. ^ As we will see below, to the best 
of our knowledge either the former or the latter require- 
ments have been enforced in the derivation of numerical 
schemes, but never both of them. The algorithms intro- 
duced here will improve both the algorithms of |6| and 
of [H. 



II. BRIEF REVIEW OF THE BENCHMARK 
ALGORITHMS AND SOME DEFINITIONS 

To assess how well each algorithm is performing, 
we start from the knowledge that for V(x) which are 
bounded from below, Eq. pleads to an equilibrium dis- 
tribution P(x,v) for the variables x and v of the form 



P(x,v) = iVexp{- (v 2 /2 + V(x)) /T). 



(2) 



where TV is a normalization constant. We are going to 
compare the exact theoretical equilibrium distribution 
to the equilibrium distribution obtained from the sim- 
ulations. It is possible, in principle, to check theoreti- 
cally which is the equilibrium distribution which is ex- 
pected integrating using a given numerical scheme, fol- 
lowing 17]: suppose we have a numerical scheme of the 
form 

Xi(t + h) = Xi(t) + hWi(xi,£) 
then the probability distribution of Xi satisfies 

00 d d 
P(xi,t + h)-P(xi,t) = ^^ — ...—K 1 ... n P(x i ,t) 



with 



n— 1 Xi 



K x ... n = {-l) n -AW x ...W n )^ 

TV. 



(3) 



2 



where (. . .)j means averaging over the noise realizations. 
In general, for systems in detailed balance at temperature 
T, 



P(Xi,t = oo) s 



P{x l ,Oo)true X exp 




Sn/T 



where P S i m is the equilibrium distribution generated in 
the simulations, and Ptrue is the theoretical equilibrium 
distribution. Given the explicit form of the numerical 
scheme, the various Kx_.. n can be computed: applying 
Eq. and expanding the r.h.s. of Eq. in the small 
parameter h, assuming (equilibrium) that P(xi, t + h) — 
P(Xi,i) = 0, we can derive the equations satisfied by 
the Si. The first nonzero Si yields the correction to the 
true equilibrium distribution generated by the numerical 
scheme. 

Let us show how to use Eq. [3] taking one of the algo- 
rithms of [(J. This algorithm integrates Eq. fusing the 
prescription 



x(t + h) = x(t) + ahv(t) + c 2 h 2 F(x(t)) - 
v(t + h) = c v(t) + dhF(x(t)) + Tj 2 



m 



(4) 



where 
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and where r\\ and r\i are two random Gaussian variables 
with zero average and moments 



(vl) = — 2 



3 - 4e-T /l 



-2-yh 



(V: 



Til 



7 h 



-2jh\ 



This implies that this algorithm fails to reproduce the 
correct equilibrium distribution at 0(h) in the exponent. 
It is possible, for the case when V[x) = w 2 x 2 /2, to solve 
this partial differential equation, obtaining the numerical 
equilibrium distribution at lowest order in h, which reads 

P(x,v) = TVexpj- [v 2 /2 + lu 2 x 2 /2] /f\ 



with 



T = 



T 



1 



27 



i.e. the numerical equilibrium distribution is similar to 
the correct one, but with a renormalization of the tem- 
perature. In particular, this effective temperature (the 
temperature "simulated" by the algorithm) goes to zero 
in the limit 7^0. In 6] it is acknowledged that the 
algorithm does not work well in this limit, although no 
formal proof is provided. 

To overcome the problem with the case of small 7, in 
a second algorithm is proposed, which reads 

x(t + h) = x(t) + C\hv(i) + C2h 2 F(x(t)) + 771 
v(t + h) = c v(t) + (ci - c 2 )hF(x{t)) + 

c 2 hF(x(t + h)) + T)2. (5) 

Using Eq. [3]to evaluate the correction to the true equi- 
librium distribution generated by this algorithm, we find 
that the contribution Si vanishes, and we are left with 
the term In other words, this algorithm reproduces 
the correct equilibrium distribution at 0(h), but there 
are still corrections 0(h 2 ) in the exponent. The algo- 
rithm given in Eq. 0is the reference algorithm which we 
propose to improve in the next section. We will refer to 
this algorithm as "Li" (from Liquid) in the following. 



(vim) = - (1 - e-^'f 

7 

The algebra to derive the correction to the equilibrium 
distribution induced by the numerical scheme in the gen- 
eral case and for a given high order integration scheme 
can be formidable; however, for a flow like Eq. ^and for 
the scheme given by Eq. the algebra is manageable. 
Assuming that (see Eq. |3 with S = Si ) 

P(x,v) = iVexp{- [v 2 /2 + V(x) +hS(x,v)] /T) , 

plugging the scheme of Eq. 0] into Eq. [31 we have that 
S(x,v) satisfies the partial differential equation 

d 2 S(x,v) _ v dS{x,v) _ ( F(x) _ v \ dS{x,v) 
dv 2 Tj dx V T l TJ dv 



F'(x) F'(x) = 



III. QUASI SYMPLECTIC ALGORITHMS 

A symplectic algorithm is a numerical scheme which 
attempts to preserve the 2-forms dqi x dpi during the 
integration of a Hamiltonian flow. The quantity qi is a 
generalized coordinate, and pi is the corresponding con- 
jugate momentum. A nice introduction to the symplectic 
integration can be found in [iH Hfl| . Given the Hamilto- 
nian H(qi,pi), and the equations of motion 

qi = {qi,H} pi = {pi,H}, 

a symplectic integrator will in practice conserve some 
quantity H , which in general reads 

H = H + h n G(pi,qi) 

where h is the integration time step, and G is a function 
which depends on the numerical scheme used for the in- 
tegration. The problem of Hamiltonian flows in the pres- 
ence of fluctuations has been addressed also in [lU [lj , 
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whereas quasi symplectic schemes were derived in ^| 
(see also below, when various comparisons are carried 
out). A preliminary account of the material of this sec- 
tion can be found in [l6| . 

Given that we are interested here in the integration 
of Eq. ^ we start from the symplectic integration of 
Hamiltonians which are separable and quadratic in the 
velocities. There are very many different possible sym- 
plectic schemes: however, having in mind that we are 
seeking a scheme which should be used in the integration 
of a stochastic differential equation, we restrict ourselves 
to considering a scheme in the form 

q(i) = q(i - 1) + ha % p(i - 1) 



p{i) = p(i - 1) + hbiF{q{i)) 

for i between 1 and N , where g(0) = q(t = 0), q(N) — 
q(t = h), etc., and h is the integration time step in the 
simulations. The coefficients a(i) and b(i) are chosen as 
to minimize, in some sense, the quantity G{p, q). 

The lowest possible symplectic algorithm one can write 
to integrate Eq. ^ following this approach reads when 
7 = T = (this scheme is also known as "leap frog" ) 



x 

v(t + h) 
x(t + h) 



= v(t) + hF(x) 
= x + -v(t + h) 



where x is the position and v is the velocity. This scheme 
conserves the quantity H — h 3 (vFF' + u 3 F"/6) / 4, where 
H = v 2 /2 + V(x). It is then possible to reintroduce both 
the dissipation (—jv) and the noise, writing the tentative 
scheme 



x = x(t) + -v(t) 
v(t + h) = c 2 [ctv(t) + hF(x) + dtTj] 
x(t + h) = x+ ~v(t + h) 



(7) 



well behaved in the limit of 7 — > 0; it has the same accu- 
racy in computing the equilibrium distribution as of Eq. 
El but it requires only one random deviate per integra- 
tion time step (as opposed to two deviates for the scheme 
of Eq. EJ), so it will run faster. In the following, we will 
refer to the algorithm of Eq. [71 as "SLO" (Symplectic 
Low Order). 

Looking at the structure of the previous algorithm, 
we can try to derive an algorithm of higher order. In 
the derivation of Eq. [7J when we applied Eq. El given 
the number of unknown quantities we could only impose 
that the term 0{h) in the exponent disappeared. If we 
could somehow increase the number of unknown quanti- 
ties when applying Eq. El while keeping the algorithm 
simple, we might be able to make the terms 0(h 2 ) in the 
exponent disappear. We start combining two steps, each 
one of the form of Eq. UJ done with an integration time 
step h/2, 



Vi 

v{t + h) 
x(t + h) 



x(t) + -v(t) 



(-2 



Cl 



v(t) + hF(xi) + y / 'jTh(airi 1 + 02^2) 



xi + -v\ 



= c 2 



V\ + hF(x2) + v ^Th(bir]i + 62772) 



X2 + -rv(t + h) 
4 



(9) 



(0) where c\ = (1 - 7V 4 ) and c 2 = 1/(1 +7/1/4). Here, 771 
and 772 are two random Gaussian variables of standard 
deviation one and average zero. The idea is now to choose 
the coefficients a's and 6's in such a way as to annihilate 
Si, and possibly minimize S2. This is done using Eq. EJin 
Eq. El which results in a number of algebraic equations 
for cii and b{. The algebra, although straightforward, is 
cumbersome and we will simply report here the results. 
For a given 02, the following choice for the other three 
parameters will ensure that Si vanishes identically: 



bo = - 



(12 



2y/U- 12a 



where 77 is a Gaussian variable, with standard deviation 
one and average zero. We use again Eq. El and, imposing 
that the term 0(h) in the exponent (i.e. the term hSi) 
vanishes, we find that the unknown arbitrary quantities 
ci, C2 and di read 



Cl 
C2 

d x 



1 



7 h 

~2 
1 



1 + 7V 2 



(8) 



Although we will carry out more extensive comparisons 
further down, let us briefly compare this scheme to the 
scheme of Eq. El The present scheme is by construction 



2%2a 2 2 + 2402^14- 16c 



42 



«i 



426i(-7V2 + 6V2al + 24a 2y /7 - 6a| 



V3(-14 



588a2) 



As function of 02, we can now write the equations sat- 
isfied by S2: we find that S2 vanishes for a particular 
choice of the parameter 02- Summarizing the numerics, 
the set of a's and 6's which simultaneously makes 5i and 
S2 vanish are: 



«i 



4.0691860043307065. 
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a 2 = -0.1533230407019893. 
h = 0.3044913128854065.. 
b 2 = -1.0363164126095790. 



(10) 



The conclusion is that the algorithm given by Eqs. El 
and 1101 is symplectic in the limit 7-^0 (conserving the 
quantity H — h 3 (vF F' + v 3 F" /6) /16) , whereas for afinite 
7 it reproduces the correct equilibrium distribution with 
an error 0(h 3 ) in the exponent. We will refer to this 
algorithm as to "SHO" (Symplectic High Order). 

In the following, we will use also the Heun algorithm to 
carry out the various comparisons. To make this paper 
as self contained as possible, we recall here that the Heun 
algorithm for a system like the one in Eq. \I]is given by 
the prescription: 

X\ = x(t) + hv(t) 

vi = v(t) - hyv(t) + hF{x(t)) + v * / 2jThr) 

%2 = x(t) + hvi 

v 2 = v(t) - hjvi + hF(xi) + ^2-yThr) 

u , h\ Xl +X2 n , h\ v i + v 2 
x (t + h) = — v(t + h) = — - — 

The Heun algorithm does not make use of the (quasi) 
symplectic nature of the flow: we expect that it will not 
fare too well in the limit of small 7. We recall that it 
is known [I| that the equilibrium distribution generated 
by the Heun algorithm is accurate up to o(h 2 ) in the 
exponent. We will refer to the Heun scheme as to "He". 

We will also compare our algorithms to the quasi sym- 
plectic algorithms of |l3j |: we should mention here that 
really the latter are weak integration schemes (for a def- 
inition of weak and strong integration schemes, see 0), 
hence they are bound to give worse results than the other 
schemes when, as we do, average quantities are consid- 
ered. The two alg orithms considered integrate with the 
prescriptions ((FJ should be consulted for more details): 



having said this, if we used Eq. [3]to asses these two algo- 
rithms, we would find that they both have a correction to 
the equilibrium distribution O(h) in the exponent: these 
will reflect in the numerical experiments, as we will com- 
ment below. 

Finally, it should be noted that given the structure of 
the Hamiltonian in the limit T — ► and 7 — > 0, which is 
H = v 2 /2 + V(x), the equilibrium distribution for Eq. 
can also be written as 

P(x,v) cx exp — H/T. 

At first sight, it may appear that the request of a sym- 
plectic integration scheme is redundant, once we made 
sure that the "correct" equilibrium distribution is gen- 
erated in the numerical integration. This is not right: 
the limit T — > is singular, hence a symplectic form for 
the numerical scheme can (and should) be imposed as an 
additional condition. 
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FIG. 1: Result of simulations for different integration 
schemes, as function of the integration time step h, for the 
system in Eq. 1141 Various moments are considered (see text 
for details), taking T = 0.1 and 7 = 1 



MTI: 



x(t + h) = x(t) + hv(t + h) 

v(t + h) = v(t) - hV'{x(t + h)) - hjv(t + h) 

+ ^2Th 1V (11) 

where v(t + h) and x(t + h) should be found recur- 
sively. 



MT2: 



v(t + h) = (1 - jh){v(t) - hV'{x(t)) + ^/2Thy^)) 
x(t + h) =x(t)+h(v{t)-hV'(x(t))) (12) 



OOHc 
□ □ SLO 
<> O SI 10 
A A Li 
<1 < MTI 
V V MT2 



< < < < 



7 V □ □ 

i X & & &-\ 



0.095 

1 

0.033 



<! O < A A' 



5*- 



□ § a v 



The random variables rj take the values ±1. These vari- 
ables are faster to generate than a Gaussian variable, 
hence these algorithms will run faster, allowing for a 
smaller integration timestep to compensate for less accu- 
racy when averaged quantities are considered. However, 



FIG. 2: Result of simulations for different integration 
schemes, as function of the integration time step h, for the 
system in Eq. 1141 Various moments are considered (see text 
for details), taking T = 0.1 and 7 = 10" 2 
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IV. NUMERICAL EXPERIMENTS 

We compare now the results of applying the algorithms 
described previously to the integration of two prototype 
stochastic differential equations. Let us first consider how 
the different algorithms reproduce the equilibrium prop- 
erties: to this end, we integrate the equations of motion, 
and compute some equilibrium momenta, which are then 
compared to the theoretical ones. 

The first system studied is given by 



v = — jv — V'(x) + \p2r{Tr] 
V(x) = x i /A-x 2 /2 



(13) 



For this system we fixed the noise intensity to T = 0.1, 
and carried out the numerical integration for two differ- 
ent values of the damping coefficient 7 and for the dif- 
ferent integration schemes. The results are summarized 
in Fig. [2 (for 7=1) and in FigH (for 7 = 1CT 2 ). The 
quantity (H) is defined as (H) = (v 2 /2 + V(x)). In all 
figures, the results of the digital simulations are shown 
by symbols with a gray straight line as guide to the eye; 
the bold dashed line is the expected (theoretical) value 
for the quantity considered. For the number of averages 
considered, the statistical error is much smaller than (or- 
der of) the symbols for the larger (smaller) damping. 

Let us comment the results. It is evident that the Heun 
method (He in figures) is not very appropriate for the 
smaller damping considered (Fig.|2J). Even for the larger 
damping (Fig. the Heun algorithm is typically out- 
performed by the Symplectic Low Order scheme (SLO, 
Eq. UJ; note that the SLO is fairly faster than He, given 
that it requires only one evaluation of the deterministic 
force for each integration time step. 

The algorithms MT1 and MT2, as expected, do not 
work well for the larger damping considered, and become 
more accurate as the damping is reduced: it should be 
noted that for this case, MT2 seems to be more accurate 
than MT1 for a given integration time step: considering 
that MT2 is much faster than MT1, the conclusion seems 
to be that MT2 ought to be preferred, between these 
two schemes. Note also that the error on the moments 
for these two schemes seems to grow linearly with the 
integration time step h, which is related to the 0(h) error 
in the exponent which was mentioned in the previous 
section. 

The Li algorithm is not particularly accurate when the 
x 2 moments are considered: for both values of the damp- 
ing, SLO gives more accurate results for this moment, in 
the whole h region. The Symplectic High Order (SHO) 
is the algorithm which gives the most accurate results for 
the x 2 moment, and results comparable or better than to 
the one obtained with Li for the v 2 and v A moments. It is 
only when (H) is considered, and for the larger damping, 
that Li seems to be more accurate than SHO. However, 
care is necessary in drawing conclusions from (H) : look- 
ing for instance at Fig. ^ we note that Li underestimates 



both v 2 and x 2 : recalling the structure of the potential 
V(x) — —x 2 /2 + x 4 / , it is clear that these two under- 
estimates tend to cancel out, leading to a (H) closer to 
the theoretical one, but only by virtue of a coincidental 
cancellation. 
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FIG. 3: Result of simulations for different integration 
schemes, as function of the integration time step h, for the 
system in Eq. 1151 Various moments are considered (see text 
for details), taking T = 0.1 and 7 = 1 
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FIG. 4: Result of simulations for different integration 
schemes, as function of the integration time step h, for the 
system in Eq. 1151 Various moments are considered (see text 
for details), taking T = 0.1 and 7 = 10 -2 

The second system studied is similar to the first one: 



x = v 

i) = — "fv — V'(x) + \j2rfTri 
V(x) = x 4 /4 + x 2 /2 



(14) 



the only difference with the system of Eq. [21 being that 
now the potential is monostable. 

The result of the computer experiments are summa- 
rized in Figs. [SJand 0| The comments parallel the com- 
ments we already made for the system of Eq. 1141 Heun 
(He) is the least accurate scheme for small damping, al- 
though the error on the moments is quadratic on h (a 
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signature of an 0(h 2 ) error in the equilibrium distribu- 
tion). MT1 and MT2 perform better at smaller damping, 
with an error on the moments which is roughly linear in 
the integration time step. SLO does better than both He 
and MT1, MT2, and for both damping considered, being 
as fast (if not faster) than both schemes. When the x 2 is 
considered, Li appears to perform worse than SLO. SHO 
outperforms Li: only when H is considered, Li seems to 
be more accurate than SHO, but again only by virtue of 
a cancellation (again, between (a: 2 ) and (i> 2 )). 
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FIG. 5: Result of simulations for different integration 
schemes, as function of the integration time step h, for the 
system in Eq. 1141 The Mean First Passage Time between the 
minima is considered (see text for details), taking T — 0.1 
and 7 = 1 
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summarized in Figs. |3]and|Bl Given the fairly large value 
of the noise intensity (T = 0.1), there is no theory avail- 
able to compute an "exact" MFPT for comparison with 
the simulated MFPTs. To have some reference, we took 
as reference value the average of the MFPT obtained via 
the algorithms Li and SHO for the three smallest h values 
used in the simulations, and drew the dashed line at this 
value. The expected statistical error, due to the finite- 
ness of the sample used, is of the order of the symbol 
size. For the larger damping considered (Fig. 7 = 1), 
He and SLO perform in a similar way. MT1 and MT2 
give unreasonable values for the MFPT (only one point 
for MT1 is actually on the graph, for the smallest h con- 
sidered: all other points for both algorithms are outside 
the MFPT range considered). Li and SHO perform in a 
similar way, giving results closer to the "correct" MFPT 
throughout the h range considered, with a slightly better 
agreement shown by SHO for larger h's. The situation is 
more interesting when a smaller damping is considered 
(Fig. H 7 = 10~ 2 ). While showing an error which grows 
only quadratically in h, clearly He is the algorithm which 
performs worse. MT1 and MT2 now give more reason- 
able results, and they yield MFPT comparable to the 
ones obtained using SLO. SHO outperforms Li, giving 
MFPT closer to the "exact" ones over the whole h range 
considered. Li on the other hand seems to give results 
which are roughly equivalent to the ones obtained using 
MT1, MT2 or SLO. 

We would like to conclude noting that the numerical 
experiments were done stretching the algorithms into pa- 
rameter regions which are somehow extreme: the typical 
time scale for the potential considered is around 0.5 (the 
oscillation frequency around the minima for Eq. 114(1 or 
around 1 (the larger 7 considered, and the oscillation fre- 
quency for the potential of Eq. I15fl , and yet an algorithm 
like SHO is able to integrate up to integration time steps 
h order of 0.3, with corrections to the moments or to the 
MFPTs which are smaller than, or at worse order of 1%: 
in our opinion, these are remarkable results, particularly 
when the flatness of the MFPT computed with SHO in 
Fig. El is considered. 



CONCLUSIONS 



FIG. 6: Result of simulations for different integration 
schemes, as function of the integration time step h, for the 
system in Eq. 1141 The Mean First Passage Time between the 
minima is considered (see text for details), taking T — 0.1 
and 7 = 10" 2 

We turn now to some numerical experiments to as- 
sess how the various algorithms perform when dynamical 
quantities are considered. Using the system of Eq. [21 
we computed the Mean First Passage Times (MFPT) to 
go from one of the minima to the other one (the minima 
are located at x = ±1): the results of the simulations are 



We introduced two algorithms for the numerical inte- 
gration of the equations of motion of a Brownian walker. 
The features of these algorithms are that they become 
symplectic when the damping on the Brownian walker is 
taken to be zero, and give the correct equilibrium dis- 
tribution to some higher order in the integration time 
steps for a finite damping and temperature. This, in 
turn, leads to more accuracy when dynamical quantities 
are considered (like the MFPT). Possible applications of 
these algorithms, beside the mentioned generic integra- 
tion of the dynamics in the liquid state 0, are in the 
integration of the dynamical equations of ions moving in 
and around ionic channels |l8[ : here the speed up pro- 
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vided by algorithms which are stable for fairly large time steps may help in improving current simulations. 
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